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ABSTRACT 

We model the long-term evolution of the Hilda collisional family located in the 3/2 
mean-motion resonance with Jupiter. Its eccentricity distribution evolves mostly due 
to the Yarkovsky/YORP effect and assuming that: (i) impact disruption was isotropic, 
and (ii) albedo distribution of small asteroids is the same as for large ones, we can 
estimate the age of the Hilda family to be 4l° Gyr. We also calculate collisional activity 
in the J3/2 region. Our results indicate that current collisional rates are very low for 
a 200 km parent body such that the number of expected events over Gyrs is much 
smaller than one. 

The large age and the low probability of the collisional disruption lead us to the 
conclusion that the Hilda family might have been created during the Late Heavy 
Bombardment when the collisions were much more frequent. The Hilda family may 
thus serve as a test of orbital behavior of planets during the LHB. We tested the 
influence of the giant-planet migration on the distribution of the family members. 
The scenarios that are consistent with the observed Hilda family are those with fast 
migration time scales ~ 0.3 Myr to 3 Myr, because longer time scales produce a family 
that is depleted and too much spread in eccentricity. Moreover, there is an indication 
that Jupiter and Saturn were no longer in a compact configuration (with period ratio 
Ps/Pj > 2.09) at the time when the Hilda family was created. 
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1 INTRODUCTION 

There are many independent lines of evidence that the or- 
bits of planets of the Solar System were not the same all the 
time, but that they have changed substantially over billions 
of years. The major arguments are based on the observed 
orbital distribution of Kuiper belt objects (Malhotra et al. 
1995, Levison et al. 2008) or small but non- negligible eccen- 
tricities and inclinations of the giant planets (Tsiganis et al. 
2005). Observations of Jupiter's Trojans (Morbidelli et al. 
2005), main-belt asteroids (Minton & Malhotra 2009, Mor- 
bidelli et al. 2010), the amplitudes of secular oscillations of 
the planetary orbits (Morbidelli et al. 2009, Brasser et al. 
2009), or the existence of irregular moons (Nesvorny et al. 
2007) provide important constraints for planetary migration 
scenarios. 

Asteroids are a fundamental source of information 
about the evolution of the planetary system. Some of the 
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resonant groups, i.e., located in the major mean-motion res- 
onances with Jupiter, might also have been influenced by 
planetary migration, because their current distribution does 
not match the map of the currently stable regions. For in- 
stance, there are two stable islands denoted A and B in the 
J2/1 resonance and only the B island is populated (Nesvorny 
& Ferraz-Mello 1997). 

In this work we focus on the Hilda asteroid family in 
the 3/2 resonance with Jupiter. We exploit our ability to 
model long-term evolution of asteroid families, which is usu- 
ally dominated by the Yarkovsky effect on the orbital ele- 
ments (Bottke et al. 2001), often coupled to the YORP effect 
on the spin rate and obliquity (Vokrouhlicky et al. 2006b). 
Chaotic diffusion in eccentricity and sometimes interactions 
with weak mean-motion or secular resonances (Vokrouhlicky 
et al. 2006a) also play important roles. In case of asteroids in- 
side strong mean-motion resonances, one has to account for 
the "resonant" Yarkovsky effect, which causes a systematic 
drift in eccentricity (Broz & Vokrouhlicky 2008). This is dif- 
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ferent from usual non-resonant orbits where the Yarkovsky 
effect causes a drift in semimajor axis. 

The Hilda collisional family — a part of the so called 
Hilda group in the 3/2 mean motion resonance with Jupiter 
— was already briefly discussed by Broz & Vokrouhlicky 
(2008). However, the modelling presented in that paper was 
not very successful, since the resulting age of the family 
seemed to be too large (exceeding 4Gyr). This was an im- 
portant motivation for our current work. We think that 
we missed an important mechanism in our previous model, 
namely perturbations arising from the migration of the gi- 
ant planets and also an appropriate treatment of the YORP 
effect. Indeed, the age > 4Gyr suggests that the planetary 
migration might have played a direct role during the early 
evolution of the Hilda family. In this paper we thoroughly 
test this hypothesis. 

The paper is organised as follows: at first, we study the 
observed properties of the J3/2 resonance population in Sec- 
tion [21 Our dynamical model of the Hilda family (without 
migration first) is described in Section [31 Then we estimate 
the collisional activity in the J3/2 region in Section [4l The 
results of our simulations of the giant-planet migration are 
presented in Section [SI Finally, Section [SI is devoted to con- 
clusions. 



2 CURRENT ASTEROID POPULATION IN 
THE J3/2 RESONANCE 

Asteroids located in the 3/2 mean motion resonance with 
Jupiter have osculating semimajor axes around (3.96 ± 
0.04) AU, i.e. beyond the main asteroid belt. Contrary to 
the Kirkwood gaps (associated with J3/1, J7/3 or J2/1 res- 
onances), this resonance is populated by asteroids while its 
neighbourhood is almost empty. The Hilda collisional family 
we are going to discuss in detail is a small part of the whole 
J3/2 resonant population. 

Our identification procedure of the J3/2 resonant popu- 
lation was described in the previous paper Broz & Vokrouh- 
licky (2008). Using the AstOrb catalogue of orbits (version 
JD = 2455500.5, Oct 31st 2010) we identified 1787 num- 
bered and multi-opposition bodies with librating critical ar- 
gument 
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where p = 2, 5 = 1, A' is the mean longitude of Jupiter, A 
the mean longitude of the asteroid and vo the longitude of 
perihelion of the asteroid. 

In order to study the detailed distribution of the bodies 
librating inside the resonance we have to use pseudo-proper 
resonant elements defined as approximate surfaces of sec- 
tions (Roig et al. 2002), i.e. intersection of the trajectory 
with a plane defined by: 
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These conditions correspond to the maximum of the semi- 
major axis a over several oscillations and the minimum of 
the eccentricity e or the inclination /. We need to apply a 
digital filter to a(t') prior to using Eq. ([2]), namely filter A 
from Quinn, Tremaine & Duncan (1991), with sampling lyr 
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Figure 1. The number A'^ of the Hilda family members versus 
the selected cut-ofF velocity iicutoff ■ 



and decimation factor of 10, to suppress fast ~ 80 yr os- 
cillations, which would otherwise disturb slower ~ 280 yr 
oscillations associated with resonant librations. Finally, we 
apply an averaging of the sections a, e, 7 over 1 Myr running 
window and these averages are the pseudo-proper elements 

Qp , Bp , 7j 

the order lO^** 
much smaller than the structures we are interested in. 

The overall dynamical structure of the J3/2 resonance is 
determined by secular resonances , at high eccentricities 
Bp > 0.3 and secondary resonances at lower values of Cp < 
0.13 (according to Morbidelli & Moons 1993, Nesvorny & 
Ferraz-Mello 1997, Ferraz-Mello et al. 1998, Roig & Ferraz- 
Mello 1999). They destabilise the orbits at the borders of a 
stable island. The orbits inside the island exhibit very low 
chaotic diffusion rates, so bodies can remain there for 4 Gyr 
(without non-gravitational perturbation). 

Next we apply a hierarchical clustering method (Zap- 
pala et al. 1994) to detect significant clusters. We use 
a standard metric in the pseudo-proper element space 
(ap,ep, sin/p) 
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In the following, we do not discuss the known Schubart 
family, which was sufficiently analysed elsewhere (Broz & 
Vokrouhlicky 2008), but we focus on the family associated 
with (153) Hilda. A suitable cut-off velocity for the Hilda 
family seems to be Wcutoff = 140 m/s, because the number 
of members does not change substantially around this value 
(see Figure[T|. The number of members at this cut-off is 400. 

The resulting plots {ap,H), {ep,H) and {Ip,H) of the 
Hilda family show very interesting features (see Figure [2} . 
The distribution of semimajor axis and inclination seems 
rather uniform and almost independent of absolute mag- 
nitude H, but eccentricities of small asteroids (i.e., with 
high H) are clearly concentrated at the outskirts of the fam- 
ily and depleted in the centre. 

In order to explain the distribution of asteroids in the 
(cp, H) plane we have to recall that asteroids orbiting about 
the Sun are affected by non-gravitational forces, mostly by 
the Yarkovsky/YORP effect, i.e. the recoil force/torque due 
to anisotropic emission of thermal radiation. We consider 
the concentrations in the {ep,H) plane to be a strong in- 
dication of the ongoing Yarkovsky/YORP evolution, be- 
cause they are very similar to those observed among sev- 
eral main-belt families in the (up, H) plane and successfully 
modelled by Vokrouhlicky et al. (2006b). The difference be- 
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tween these two cases stems from the fact that main-belt 
famiHes are non-resonant and the Yarkovsky/YORP effect 
thus increases or decreases the semimajor axis (depending on 
the actual obliquity of the spin axis), while in our resonant 
case, the same perturbation results instead in a systematic 
increase or decrease of eccentricity. A detailed modelling of 
the e-distribution is postponed to Section [3.51 

The central part of the {ep,H) distribution, from e = 
0.17 to 0.23, seems rather extended. The large asteroids 
{H < 12.5 mag) are spread over this interval of eccentrici- 
ties even though their Yarkovsky drift rates must have been 
small. Only 2-4 of them are likely to be interlopers, be- 
cause there is a very low number of background asteroids 
in the surroundings of the family (see Figure [3]). We think 
this shape might actually be the result of the initial size- 
independent perturbation that the family distribution re- 
ceived by the migration of the giant planets (which we dis- 
cuss in Section [5. ip . 

Regarding the {ap,H) distribution, the largest asteroid 
(153) Hilda is offset with respect to the centre, but this is 
a natural outcome of the definition of the pseudo-proper 
elements — fragments which fall to the left of the libration 
centre are mapped to the right which creates the offset. 

The geometric albedos for Ifilda family objects are 
poorly known. There are only six measured values for the 
family members: 0.064, 0.046, 0.038, 0.089, 0.044, 0.051 
(Davis & Neese 2002). Given the low number of values and 
the possibility of selection effects we prefer to assume the 
family members have a mean value pv = 0.044, which corre- 
sponds to the whole J3/2 population. The size of the parent 
body can be then estimated to be DpB = (200 ± 20) km. We 
employ two independent methods to determine the diameter 
DpB'- (i) we sum the volumes of the observed bodies larger 
than an assumed completeness limit Dcompiete = 10 km and 
then we prolong the slope of the size-frequency distribu- 
tion down to -D = to account for unobservable bodies (see 
Broz & Vokrouhlicky 2008), which results in Dpb — 185 km; 
(ii) we also use a geometric method developed by Tanga et 
al. (1998) which gives Dpb — 210 km. A test with different 
albedo values will be described in Section [3.61 

The size-frequency distribution N{>D) vs D of the 
Ifilda family is steeper than that of background J3/2 popu- 
lation, but shallower than for usual main-belt families (Fig- 
ure |4|. Interestingly, the slope 7 — —2.4 ±0.1 of the distri- 
bution N{>D) = CD'' is close to a collisional equilibrium 
calculated by Dohnanyi (1969). 

Colour data extracted from the Sloan Digital Sky Sur- 
vey Moving Object Catalogue version 4 (Parker et al. 2008) 
confirm the Hilda family belongs to the taxonomic type C, 
because most of the spectral slopes are small. Recall that 
the whole J3/2 population exhibits a bimodal distribution 
of slopes, i.e. it contains a mixture of C- and D-type aster- 
oids. 



3 THE HILDA FAMILY MODEL WITH 
RADIATION FORCES 

To understand the long term evolution of the Hilda family, 
we construct a detailed numerical model, extending efforts 
in Broz & Vokrouhlicky (2008), which includes the following 
processes: (i) impact disruption, (ii) the Yarkovsky effect, 
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Figure 3. The J3/2 region displayed in (ap,7p) plot. A very 
prominent Schubart cluster (studied by Broz & Vokrouhlicky 
2008) is visible around sin/p = 0.05. The close surroundings of 
the Hilda family, where only a low number of bodies is present, 
is highlighted by grey rectangles. 
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Figure 4. Cumulative size distributions of the J3/2 population 
and the Hilda family. The polynomial fits of the form N{>D) = 
CD"/ are plotted as thin lines, together with the respective val- 
ues of the 7 exponent. Several main-belt families are plotted for 
comparison: Eos (with slope 7 = —2.8), Eunomia (—5.0), Hygiea 
(-3.8), Koronis (-2.8), Themis (-2.9), Tirela (-3.3), Veritas 
(-3.4) and Vesta (-5.4). 



(iii) the YORP effect, (iv) collisions and spin-axis reorienta- 
tions. We describe the individual parts of the model in the 
forthcoming subsections. 



3.1 Impact disruption 

To obtain initial conditions for the family just after the 
breakup event we need a model for the ejection velocities of 
the fragments. We use a very simple model of an isotropic 
ejection from the work of Farinella et al. (1994). The distri- 
bution of velocities " at infinity" follows the function 

-(c«+l)/2 J 



dN(v)dv = Cv{v^ +vi^y 



(4) 



with the exponent a being a free parameter, C a normali- 
sation constant and Vasc the escape velocity from the parent 
body, which is determined by its size Dpb and mean den- 
sity ppB as VcBc = \/ (2/3)7rGppB Dpb ■ The distribution is 
usually cut at a selected maximum allowed velocity Umax to 
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Figure 2. The Hilda family displayed in resonant semimajor axis ap (left), eccentricity ep (middle) and inclination sin/p (right) versus 
absolute magnitude H. The libration centre is located at a ~ 3.96 AU and all bodies are displayed to the right of it. The 'ears' in (ep, H), 
i.e., the concentration of small asteroids at the outskirts of the family and their depletion in the centre, are very prominent here. The thin 
vertical lines denote the central part of the (ep, J?) distribution discussed in the text. The family has 400 members at ?^cutofI = 140 m/s. 



prevent outliers. The actual values of all these parameters 
are given in Section 13.51 Typically, the overall distribution 
of velocities has a peak close to the escape velocity, which is 
approximately 100 m/s for a 200 km parent body. The ini- 
tial velocities \v\ of individual bodies are generated by a 
straightforward Monte-Carlo code and the orientations of 
the velocity vectors v in space are assigned randomly. 

Here, we assume the velocity of fragments is indepen- 
dent of their size, which seems reasonable with respect to 
the observed uniform distribution of the Hilda family in the 
(flp, H) and (Jp, H) planes (Figure [S]). We perform also tests 
with non-isotropic distributions in Section [3.71 

We must also select initial osculating eccentricity ei 
of the parent body, initial inclination ii, as well as true 
anomaly /imp and argument of perihelion LJimp at the time 
of impact disruption. All of these parameters determine the 
initial shape of the synthetic "Hilda" family just after the 
disruption of the parent body. Initial semimajor axis ai is not 
totally free, instead it is calculated from the initial semima- 
jor axis of Jupiter aji and the Kepler law, since the parent 
body has to be confined in the J3/2 resonance. 



3.2 Yarkovsky effect in a resonance 

The long-term evolution of asteroid orbits is mainly driven 
by the Yarkovsky thermal effect. The implementation of the 
Yarkovsky effect in the SWIFT integrator was described in 
detail in Broz (2006). Only minor modifications of the code 
were necessary to incorporate spin rate evolution, which is 
driven by the YORP effect (see Section [3. 3p . 

The thermal parameter we use are reasonable estimates 
for C/X-type bodies: Psurf = Pbuik = 1300 kg/m^ for the 
surface and bulk densities, K = 0.01 W/m/K for the surface 
thermal conductivity, C = 680J/kg for the heat capacity, 
A = 0.02 for the Bond albedo and eiR = 0.95 for the thermal 
emissivity parameter. 

We can use a standard algorithm for the calculation of 
the Yarkovsky acceleration which results in a semimajor-axis 
drift in case of non-resonant bodies. The drift in eccentric- 
ity in case of resonant bodies arises " automatically" due to 
the gravitational part of the integrator. In Figure [S] we can 
see a comparison between the expected drift Aa in semima- 
jor axis and the resulting drift Ae in eccentricity, computed 
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Figure 5. Almost linear relation between the expected drift Aa 
in semimajor axis and the simulated drift Ae in eccentricity, com- 
puted for 360 members of the Hilda family located inside the J3/2 
resonance. 



for the Hilda family (see the explanation in Appendix A of 
Broz & Vokrouhlicky 2008). The data can be approximated 
by a linear relationship, where the departures from linear- 
ity are caused mainly by interactions of drifting orbits with 
embedded weak secular or secondary resonances. 

Note that according to a standard solar model the 
young Sun was faint (Giidel 2007), i.e., its luminosity 4Gyr 
ago was 75 % of the current Lq . We can then expect a lower 
insolation and consequently weaker thermal effects acting 
on asteroids. Since we assume a constant value of Lq in our 
code the age estimated for the Hilda family (in Section [3.511 
can be 12.5 % larger. 



3.3 YORP effect 

The magnitude of the Yarkovsky drift sensitively depends on 
the orientation of the spin axis with respect to the orbital 
plane and, to a lesser extent, on the angular velocity too. 
We thus have to account for the long-term evolution of spins 
of asteroids which is controlled by torques arising from the 
emission of thermal radiation, i.e. the YORP effect. The 
implementation of the YORP effect follows Vokrouhlicky et 
al. (2006). We assume the following relations for the rate of 
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angular velocity and obliquity 
Aui 

At uj 



200, 



(5) 
(6) 



where /- and (j-functions are given by Capek & Vokrouhlicky 
(2004) for a set of 200 shapes with mean radius _Ro = 1 km, 
bulk density po ~ 2500 kg/m^, located on a circular orbit 
with semimajor axis ao = 2.5 AU. The shapes of the Hilda 
family members are not known, so we assign one of the arti- 
ficial shapes (denoted by the index i) randomly to each indi- 
vidual asteroid. We only have to scale the /- and (;-functions 
by a factor 



c = cyorp 1 — 
ao 



Pbulk 

po 



(7) 



where a, R, pbuik are semimajor axis, radius and density 
of the simulated body, and cyorp is a free scaling param- 
eter, which can account for an additional uncertainty of 
the YORP model. Because the values of /'s and g's were 
computed for only a limited set of obliquities (with a step 
Ae = 30°) we use interpolation by Hermite polynomials (Hill 
1982) of the data in Capek & Vokrouhlicky (2004) to obtain 
a smooth analytical functions for /i(e) and gi{e). 

If the angular velocity approaches a critical value 



(8) 



t^crit = ^1 -TfGpbulk , 

we assume a mass shedding event, so we keep the orienta- 
tion of the spin axis and the sense of rotation, but we reset 
the orbital period P = 27r/a; to a random value from the 
interval (2.5, 9) hours. We also change the assigned shape to 
a different one, since any change of shape may result in a 
different YORP effect. 

The differential equations (O, ((U are integrated nu- 
merically by a simple Euler integrator. The usual time step 
is At = 1000 yr. An example of the results computed by 
the spin integrator for the Hilda family is displayed in Fig- 
ure [6] The typical time scale of the spin axis evolution is 
tyorp — 500 Myr. After ~ 3 times tyorp most bodies have 
spin axes perpendicular to their orbits, what maximizes the 
Yarkovsky drift rate of eccentricity. 

3.4 Collisions and spin-axis reorientations 

In principle, collisions may directly affect the size distribu- 
tion of the synthetic "Hilda" family, but we neglect this ef- 
fect because most of the asteroids are large enough to remain 
intact. 

However, we include spin axis reorientations caused by 
collisions. We use an estimate of the time scale by Farinella 
et al. (1998) 



where B = 84.5 kyr, Pi = 5/6, P2 = 4/3, Do = 2m and 
ojQ corresponds to period P = 5 hours. These values are 
characteristic for the main belt and we use them as an upper 
limit of Trcor for the J3/2 region. Even so, the time scale is 
Treor — 3 Gyr for the smallest observable (_D ~ 5 km) bodies 
and reorientations are thus only of minor importance. Note 




Figure 6. An example of the YORP-driven evolution of obliq- 
uities (namely a ^-component of the spin axis unit vector, top 
panel) and angular velocities u) (bottom panel) for the members 
of the synthetic " Hilda" family. At the beginning, all values of u) 
were selected positive and spin axes were distributed isotropically. 
The evolution may force u) to become negative, which simply cor- 
responds to an opposite orientation of the spin axis. The scaling 
parameter was selected Cyorp = 0.33 in this run. 



that the probability of the reorientation is enhanced when 
the YORP effect drives the angular velocity uj close to zero. 



3.5 Results on Yarkovsky /YORP evolution 

We start a simulation with an impact disruption of the par- 
ent body and create 360 fragments. Subsequent evolution of 
the synthetic "Hilda" family due to the Yarkovsky/ YORP 
effect is computed up to 6 Gyr in order to estimate the time 
span needed to match the observed family even though the 
family cannot be older than ~ 4 Gyr, of course. Planets are 
started on their current orbits. A typical outcome of the 
simulation is displayed in Figure [T] 

Due to the long integration time span and large num- 
ber of bodies, we were able to compute only four simulations 
with the following values of true anomaly and YORP effi- 
ciency: 

(i) /imp = 0°, Cyorp = 0; 
(n) /imp = 180°, Cyorp = 0; 
(in) /imp = 0°, Cyorp = 1; 
(iv) /imp = 0°, Cyorp = 0.33. 

The remaining parameters were fixed: ei = 0.14, ii = 7.8°, 
i^imp = 30°, a — 3.25, Wmax ~ 300 m/s, i?pB = 93.5 km, 
PPB = 1300kg/m^ pv = 0.044. 

We are mainly concerned with the distribution of eccen- 
tricities Cp, because the observed family has a large spread 
of Bp's, while the initial synthetic family is very compact. 
For this purpose we constructed a Kolmogorov-Smirnov test 
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Figure 7. Eccentricity vs absolute magnitude plot for ttie syn- 
thetic "Hilda" family just after the impact disruption (time 
t = 0, top panel) and after 4 Gyr of evolution due to the 
Yarkovsky/YORP effect (bottom panel). There is a comparison 
with the observed Hilda family (gray dots). 



(Press et al. 1999) of the normalized cumulative distribu- 
tions iV(<e) 



max |Af(<e)syn 
0<e<l ' ^ ' ' 



N{<e)ob 



(10) 



which provides a measure of the difference between the syn- 
thetic "Hilda" family, at a given time, and the observed 
Hilda family (see Figure [8] for an example) . The results of 
the KS tests are summarized in Figure [5] (first four panels). 

There is an easy possibility to asses the sensitivity of 
results with respect to the Wmax parameter too, without the 
need to compute the simulation again. We simply select bod- 
ies fulfilling the condition v < u^ax, with w^ax ~ 200, 100 or 
50m/s, and recompute only the KS statistics for this sub- 
set. The results are plotted in Figure H] as thin lines. We can 
state values lower than Umax — 100 m/s are surely excluded. 

As a preliminary conclusion we may say that all sim- 
ulations point to a large age of the Hilda family. The gp- 
distributions are most compatible with the observed family 
for ages t — (4.0 ± 1.0) Gyr. This suggests the Hilda fam- 
ily might have experienced the giant-planet migration pe- 
riod which is dated by the Late Heavy Bombardment to 
fLHB — 3.85 Gyr (Gomes et al. 2005). The large uncertainty 
of the age stems from the fact that the runs including the 
YORP effect (cyorp 5S 0.33) tend to produce ages at a 
lower limit of the interval while the YORP-less runs (with 
Cyorp = 0) tend to the upper limit. 

3.6 Alternative hypothesis: high albedos of small 
asteroids 

We now discuss two scenarios that further reduce the mini- 
mal age of the family: (i) high albedos of small asteroids (i.e., 
larger Yarkovsky/YORP drift); (ii) strongly asymmetric ve- 
locity field after impact (like that of the Veritas family). 




0.35 



Figure 8. Normalized cumulative distributions N(<e) of ec- 
centricities for (i) the observed Hilda family, (ii) the synthetic 
"Hilda" family at time t = (just after the impact disrup- 
tion), (iii) evolved due to the Yarkovsky/YORP effect (at time 
t = 3845 Myr) . In this figure we show the best fit for the sim- 
ulation with parameters /imp = 0°, cyorp = 0.33. Note the 
'bended' shape of the observed distribution corresponds to the 
'ears' on the {ep,H) plot (Figure [2J. There is no perturbation by 
planetary migration in this particular case. 



Albedo is the most important unknowir parameter, 
which can affect results on the Yarkovsky/YORP evolu- 
tion. Fernandez et al. (2009) measured albedos of small 
Trojan asteroids and found a systematically larger values 
than for large Trojans. If we assume the J3/2 asteroids be- 
have similarly to Trojans, we may try a simulation with an 
rather high value of geometric albedo pv ~ 0.089 (com- 
pared to previous pv = 0.044). Moreover, we decrease den- 
sity pbuik = 1200 kg/m^, increase maximum velocity of 
fragments Umax = 500 m/s (though the velocity distribu- 
tion is still determined by Eq. |4]) and select true anomaly 
/imp ~ 90° to maximise the spread of Bp's. 

The KS test is included in Figure |9l panel (e). The most 
probable age is (2.3 ± 0.5) Gyr in this case. However, we do 
not think that the size-dependent albedo is very plausible 
because both large and small family members should origi- 
nate from the same parent body and their albedos, at least 
just after the disruption, should be similar. Nevertheless, the 
albedos may change to a certain degree due to space weath- 
ering processes (Nesvorny et al. 2005) . Unfortunately, we do 
not have enough data for small asteroids to assess a possible 
albedo difference between large and small family members. 



3.7 Alternative hypothesis: strongly asymmetric 
velocity field 

Another possibility to reduce estimate of the family age 
is that the original velocity was highly anisotropic. A well 
known example from the main belt is the Veritas family. 
Let us assume the anisotropy is of the order of Veritas, i.e., 
approximately 4 times larger in one direction. Note that Ver- 
itas is a young family and can be modelled precisely enough 
to compensate for chaotic diffusion in resonances (Nesvorny 
et al. 2003, Tsiganis et al. 2007). This family is character- 
istic by a large spread in inclinations, which corresponds 
to large out-of-plane components of velocities. In case of 
the Hilda family we multiply by 4 the radial components 
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Figure 9. Kolmogorov-Smirnov tests of the synthetic "Hilda" family: (a) no migration, only initial disruption (at anomaly /imp = 0°, 
roi = 30°) and subsequent Yarkovsky evolution; (b) the case with /imp = 180°; (c) including the YORP effect; (d) YORP with efficiency 
factor Cyorp = 0.33; (e) high albedo values (i.e., small bodies); (f) strongly asymmetric velocity field. The horizontal line denotes the 
distance Dks = 0.165 for which the probability p(>£'ks) that the two eccentricity distributions differ by this amount equals to 0.01. 



of initial velocities to maximise the dispersion of eccentric- 
ities, assuming the most favourable geometry of disruption 

(/imp =90°). 



The fit in Figure [9l panel (f) is seemingly better at 
the beginning of the simulation, but bodies on unstable 
orbits are quickly eliminated and the fit gets much worse 
at t ~ 500 Myr. We can see that the synthetic "Hilda" 
family is similar to the observed Hilda family quite early 
(at t ~ 2.5 Gyr), however the best fit is at later times 
[t ~ 3.5 Gyr), so there is no significant benefit compared 
to isotropic velocity-distribution cases. 



4 DISRUPTION RATES IN THE J3/2 
POPULATION 

4.1 Present collisional activity 

The results presented above show that the Hilda family is 
old. However, the uncertainty of the age is too large to con- 
clude whether the family formed during the LHB period. An 
alternative constraint is the collisional lifetime of the parent 
body. If the probability that the parent body broke in the 
last 4 Gyr in the current collisional environment is negligi- 
ble, this would argue that the family broke during the LHB 
when the collisional bombardment was much more severe. 
Thus, here we estimate the collisional lifetime of the parent 
body. 
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In our case, the target (parent body) has diameter 
ftargct = 200 km, mean impact velocity Vi — 4.8km/s 
(Dahlgren 1998), probable strength Qo = 4x 10^ J/kg (Benz 
& Asphaug 1999) and thus the necessary impactor size (Bot- 
tke et al. 2005) is 

1/3 



disrupt — {2Qiy/Vi 



Dt 



argct 



65 km . 



(11) 




100 



The population of ^65 km projectiles is dominated by main- 
belt bodies: riprojcct — 160, according to Bottke et al. (2006), 
and we have only one 200 km target in the J3/2 region, so 
ntarget = 1. The intrinsic coUisional probability for Hilda vs Figure 10. Mean intrinsic coUisional probability P; and mean 



40 60 80 

f/Myr 



main belt collisions is Pi = 6.2 x 10 '^^ km ^ yr ^ (Dahlgren 
1998) and the corresponding frequency of disruptions is 



fdl 



srupt 



Pi 



-^target 



- Tlprojcct ^target 



10 yr 



(12) 



Thus, over the age of the Solar System Tss — 4 Gyr (after 
LHB), we expect a very low number of such events 

^events — 

^SSydisrupt 

~ 0.004. 

The value of strength Ql, used above corresponds to 
strong targets. Though there is a theoretical possibility that 
the Hilda parent body was weaker, it does not seem to us 
likely, because the Hilda family is of the C taxonomic type. 
Thus, it is rather similar to (presumably stronger) main belt 
asteroids, than to (likely weaker) D-type objects. Anyway, 
even if we use an order of magnitude lower strength inferred 
for weak ice, Q|j ~ 4 x 10* J/kg (see Leinhardt & Stew- 
art 2009, Bottke et al. 2010), we obtain ddisrupt — 30km, 
Wprojcct — 360 and rievents — 0.009, so the conclusion about 
the low number of expected families remains essentially the 
same. 



4.2 The Late Heavy Bombardment 

We now compute the probability that the parent body broke 
during the LHB. We can think of two projectile popula- 
tions: (i) transient decaying cometary disk; (ii) D-type as- 
teroids captured in the J3/2. Models like that of Levison 
et al. (2009) suggest the decay time scale of the cometary 
bombardment is of the order 10 to 100 Myr and the flux of 
impactors integrated over this time span might have been 
100 times larger than today. Higher mean coUisional veloci- 
ties, due to projectiles on high-e and high-i orbits, are also 
favourable. 

In order to estimate coUisional activity we use a self- 
consistent model of the cometary disk from Vokrouhlicky, 
Nesvorny & Levison (2008). Their N-body simulations in- 
cluded four giant planets and 27,000 massive particles with 
a total mass Mdisk ~ 35 M^. The orbital evolution was prop- 
agated by the SyMBA integrator for 100 Myr. Using the out- 
put of these simulations, we calculate the mean intrinsic col- 
lisional probabilities Pi{t) between the cometary-disk popu- 
lation (at given time t) and the current J3/2 population. We 
use an algorithm described in Bottke et al. (1994) for this 
purpose. Typically, the P reaches 2 to 3 x 10~^^ km~'^yr~^ 
and the corresponding mean impact velocities are Vimp = 
7 to lOkm/s (see Figure [TU)) . 

The necessary impactor size is slightly smaller than be- 
fore, ddisrupt ~ 40 to 50 km due to larger Vimp- To estimate 
the number of such projectiles we assume that the cometary 
disk had a size distribution described by a broken power-law 
with differential slopes qi — 5.0 for D > Do, §2 = 2.5 ± 0.5 



impact velocity Vimp versus time for one of the disk simulations 
from Vokrouhlicky, Nesvorny & Levison (2008). 



for D < Do, where the diameter corresponding to the change 
of slopes is Do = 50 to 70 km. We then use the following ex- 
pressions to calculate the number of bodies larger than the 
given threshold (Vokrouhlicky, Nesvorny & Levison 2008) 



Di = Do 



N{>D) 



where Mo 



(gi -4)(4-g2) Md 



(gi - l)(gi -g2) Mq 



(13) 



92-1 

_ qi - 92 
52 — 1 

l3 



Do 



qi-l 



Do 



(14) 



fpL>5 and p = 1300kg/m'='. The result of 
this calculation is A'^(>ddisrupt ) = 0.3 to 1.7 x 10^ The ac- 
tual number of bodies in the simulation (27,000) changes in 
course of time and it was scaled such that initially it was 
equal to Ai'(>cidisrupt). The resulting number of events is 



-'-'tarect 



''events — 



^target / Pi{t) Up 

roject 

[t) dt 



0.05 to 0.2, 



(15) 



which is 10 to 50 times larger than the number found in 
Section |4T] 

Regarding the captured D-type asteroids, they were 
probably not so numerous and their impact velocities were 
lower but their coUisional probabilities were larger and the 
population might have had substantially longer timescale of 
decay (Levison et al. 2009). Using the following reasonable 
values: Vi = 4.0km/s, ddisrupt = 70 km, Tiprojoct = 5000, 
P = 2.3 X 10"^*km"^yr"^ Tlhb ^ 1 Gyr, we obtain the 
number of events ~ 0.1 which is again 25 times larger than 
the number from Section [4. II 

We conclude the Hilda family was likely created dur- 
ing the Late Heavy Bombardment when the collisions were 
much more frequent than in the current coUisional environ- 
ment. We must now test whether the structure of the family 
is consistent with the giant-planet migration, since it is con- 
nected with the LHB. 



5 PLANETARY MIGRATION 

At the LHB-time the planetary migration was most proba- 
bly caused by a presence of a massive cometary disk. Instead 
of a full N-body model we use a simpler analytic migration, 
with an artificial dissipation applied to the planets. This is 
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the only viable possibility in our case, because we need to 
test not only a large number of various migration scenar- 
ios but also various initial configurations of the synthetic 
"Hilda" family. 

For this purpose we use a modified version of the 
symplectic SWfFT-RMVSS integrator (Levison & Duncan 
1994). We account for four giant planets and include the fol- 
lowing dissipation term applied to the planets in every time 
step 

^ ^\ Av At ( t-to W 

V = V 1 H exp I I , (16) 

where v denotes a velocity vector of a given planet, v the ab- 
solute value of velocity. At the time step, Tmig the selected 
migration time scale, Av = -y/ GM/ai — \/ GM/at the re- 
quired total change of velocity (i.e., the difference of mean 
velocities between the initial and the final orbit), t the time 
and to some reference time. If there are no other perturba- 
tions than (|16p . the semimajor axis of the planet changes 
smoothly (exponentially) from the initial value ai to the fi- 
nal a{. We use time step At = 36.525 days and the total 
time span of the integration is usually equal to 3rmig when 
planetary orbits practically stop to migrate. 

We would like to resemble evolution of planetary orbits 
similar to the Nice model so it is necessary to use an ec- 
centricity damping formula, which simulates the effects of 
dynamical friction (Morbidelli et al. 2010). This enables us 
to model a decrease of eccentricities of the giant planets to 
relatively low final values. The amount of eccentricity damp- 
ing is characterised by a parameter edamp- 

Because inclinations of the planets are not very impor- 
tant for what concerns the perturbation of minor bodies (the 
structure of resonances is mainly determined by planetary 
eccentricities) , we usually start the planets with current val- 
ues of inclinations. 

We admit the analytic migration is only a crude ap- 
proximation of the real evolution, but we can use it as a 
first check to see which kinds of migration scenarios are al- 
lowed and which are not with respect to the existence and 
structure of the Hilda family. 

As a summary we present a list of free and fixed (as- 
sumed) parameters of our model in Tables [T] and O Accord- 
ing to our numerical tests the initial configuration of Uranus 
and Neptune is not very important, as these planets do not 
produce significant direct perturbations on asteroids located 
in the J3/2 resonance. We thus do not list the initial semi- 
major axes and eccentricities of Uranus and Neptune among 
our free parameters thought we include these planets in our 
simulations. 

The problem is we cannot tune all 17 parameters to- 
gether, since the 17-dimensional space is enormous. We thus 
first select a reasonable set of impact parameters for the 
family (8.-17. in Table [l|, keep them fixed, and experiment 
with various values of migration parameters (1.-7.) We test 
roughly 10^ migration scenarios. Then, in the second step, 
we vary impact parameters for a single (successful) migra- 
tion scenario and check the sensitivity of results. 

5.1 Results on planetary migration 

In the first test we compute an evolution of the syn- 
thetic "Hilda" family during planetary migration phase 
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Table 1. Free parameters of our Hilda family model. 



no. 


parameter 


description 


1. 


aji 


initial semimajor axis of Jupiter 


2. 


asi 


Saturn 


3. 


eji 


initial eccentricity of Jupiter 


4. 


esi 


Saturn 


5. 


fmig 


migration time scale 


6. 


^damp J 


eccentricity damping for Jupiter 


7. 


^dampS 


Saturn 


8. 


Ei 


initial eccentricity oi the parent body 


9. 


2i 


initial inclination 


10. 


/imp 


true anomaly at the impact disruption 


11. 




argument of perihelion 


12. 


a 


slope of the velocity distribution 


13. 


^max 


maximum velocity of fragments 


14. 




radius of parent body 


15. 


Ppb 


bulk density 


16. 


Pv 


geometric albedo of fragments 


17. 


cyorp 


efficiency of the YORP effect 



Table 2. Fixed (assumed) parameters of the Hilda family model. 
There is also a number of less important parameters, like the 
thermal ones (psurf, K, C, A, ei^) or collisional (S, /3i, /32). 



no. 


parameter 


description 


18. 


o-jf 


final semimajor axis of Jupiter 


19. 




Saturn 


20. 


N(<H) 


(observed) absolute magnitude distribution 



for the following parameter space (these are not inter- 
vals but lists of values): aji = (5.2806 and 5.2027) AU, 
asi = (8.6250,8.8250,9.3000) AU, ej; = (0.065,0.045), esi = 
(0.08,0.05), Tmig = (0.3, 3, 30, 300) Myr, edampj = 10"", 

CdampS 

= 10 El The values of a,ii and asi correspond to 
period ratios Ps/Pj from 2.09 to 2.39 (the current value is 
2.49), i.e. the giant planets are placed already beyond the 
2:1 resonance, since the 2:1 resonance crossing would de- 
stroy the Hilda family (Broz & Vokrouhlicky 2008) . Impact 
parameters were fixed except /imp: Si = 0.14, ii = 7.8°, 
/imp = (0°,180°), ojimp = 30°, a = 3.25, Umax = 300m/s, 
RPB = 93.5 km, ppb = 1300 kg/m^ 

The synthetic "Hilda" family has 360 bodies in case of 
short simulations (rmig = 0.3 or 3 Myr). In case of longer 
simulations we create 60 bodies only. Their absolute magni- 
tudes (sizes) were thus selected randomly from 360 observed 
values. This is a minimum number of bodies necessary to 
compare the distributions of eccentricities. We performed 
tests with larger numbers of bodies and the differences do 
not seem significant. 

A comparison of the final orbits of the planets with 
current planetary orbits shows we have to exclude some mi- 
gration simulations (mostly those with Uranus and Neptune 
on compact orbits). One of the reasons for unsuccessful sce- 
narios is that a compact configuration of planets is inher- 

^ In order to increase the statistics we ran simulations multiple 
times with different initial conditions for Uranus and Neptune: 
aui = (18.4479, 12.3170) AU, aw; = (28.0691, 17.9882) AU, eui = 
(0.06,0.04), CNi = (0.02,0.01). 
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Figure 11. A usual evolution of the synthetic "Hilda" family in 
the pseudo-proper semimajor axis vs eccentricity plot. The initial 
{t = Myr) and final stages (t = 100 Myr) are plotted. The mi- 
gration time scale was Tmig = 30 Myr in this particular example.. 
We selected this longer time scale because secular frequencies can 
be then computed more precisely (see Figure [12]|. The arrow in- 
dicates a total change of the position of the J3/2 resonance due 
to migration of Jupiter. 



ently unstable. If the migration time scale is too large or 
the eccentricity damping too low, it may result in a violent 
instability, close encounters between planets and eventually 
an unrealistic final configuration. 

The change in the structure of the synthetic "Hilda" 
family due to migration can be seen in Figure [TT] The fam- 
ily is shifted in semimajor axis, because it moves together 
with the resonance with migrating Jupiter. Moreover, the 
eccentricities are dispersed while the inclinations are barely 
affected. 

We identified that the eccentricity distribution is modi- 
fied when secondary resonances occur between the libration 
frequency /j3/2 of an asteroid in the J3/2 resonance and the 
frequency /ij_2S of the critical argument of Jupiter-Saturn 
1:2 resonance (see Kortenkamp et al. 2004 or Morbidelli et 
al. 2005 for case of Trojans) 



nfj3/2 = /lJ-2S : 



(17) 



where n is a small integer number, n = 2, 3 or 4 in our 
case0 We can see the evolution of resonant semimajor axes 
and the corresponding dominant frequencies, computed by 
means of periodogram, in Figure [121 

Because the resonances are localised — they act only 
at particular values of semimajor axes of planets — it is 
not necessary to have a dense grid in aji, asi parameters to 
study the dependence of the synthetic "Hilda" family shape 
on aji, asi- Essentially, there are only three situations, when 
the Hilda family is strongly perturbed, otherwise the spread 
in e does not change much in course of time. 

A very simple test, which allows us to quickly select 
allowed migration scenarios, is the number of remaining 
"Hilda" family members. We may assume the depletion by 
dynamical effects was probably low (say 50% at most), 
otherwise we would obtain much larger parent body than 
D ~ 200 km, which has much lower probability of col- 
lisional disruption. The fractions of the remaining bodies 
MeftZ-^Vinitiai vcrsus initial conditions for planets are dis- 
played in Figure [T^ 



^ We also looked for secondary resonances connected with the 4:9, 
3:7 and 2:5 Jupiter-Saturn resonances, but we found no significant 
effects. 
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Figure 12. Top panel: the frequency /ij_2S of th^ Jupiter- 
Saturn 1:2 mean motion critical argument (thick gray curve) vs 
time t. The frequency changes due to the migration of planets 
with the time scale r^ig = 30 Myr. We also computed dominant 
frequencies of librations in the J3/2 resonance for three se- 

lected members of the synthetic Hilda family (black curves). We 
do not plot the frequency itself but a selected multiple of it n/j3 /2 ■ 
Captures in the secondary resonances of type n/,j3/2 = /ij— 28 
are then clearly visible when the frequencies are equal. For the 
test particle number 1 it occurs between 4 and 10 Myr, particle 2 
was captured from 21 to 32 Myr and particle 3 from 54 Myr till the 
end of the simulation. Bottom panel: the corresponding changes 
of the pseudo-proper semimajor axes ap vs time t due to the sec- 
ondary resonances. The three test particles from the top panel 
are shown (black curves) together with the remaining members 
of synthetic "Hilda" family (gray curves). Note that some parti- 
cles may be pushed to the border of the stable libration zone and 
then escape from the J3/2 resonance. 



Low number of remaining bodies A^ioft indicates that 
perturbations acting on the synthetic family were too strong. 
It means either the family had to be formed later (when 
fewer and weaker secondary resonances are encountered) to 
match the observed family or this migration scenario is not 
allowed. The same applies to the dispersion of e-distribution 
(see below) : if it is too large compared to the observed Hilda 
family, the synthetic "Hilda" had to be formed later or the 
scenario is not allowed. Our results indicate that: 

(i) a faster migration time scale Tmig — 0.3 Myr to 30 Myr 
is preferred over slower time scales; 

(ii) Jupiter and Saturn were not in the most compact 
configuration (aji = 5.2806 AU, asi = 8.6250 AU) at the 
time when the "Hilda" family was created. 
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Figure 13. The number of simulations A'^ versus the fraction of remaining bodies A^ieft /A^initial from the synthetic "Hilda" family. The 
histograms arc plotted for four different time scales of migration r^ig and six different initial configurations of Jupiter and Saturn (aj; , 
agi; we indicate period ratios Psi/Pji instead of semimajor axes here). The ranges of remaining free parameters are mentioned in the 
text. We only plot successful migration scenarios with Aupianota ^; 2000 m/s, where Avpianots = 'Yli ^""i ^ sum of velocity differences Sv 
(defined similarly as in the HCM metric, Eq. ^ between the final simulated orbit of the j-th planet and the currently observed one. This 
way we join differences in orbital elements a, e, I into a single quantity which has the dimension of a velocity. 
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Figure 14. Eccentricity dispersions of the synthetic "Hilda" families at the end of the giant-planet migration for various initial conditions 
of the impact disruption: initial eccentricity e\, inclination ii, true anomaly /imp, argument of perihelion tJimp, exponent a, maximum 
velocity Vmax, radius of the parent body RpB and its bulk density ppB- The values of remaining parameters related to migration are 
mentioned in the text. Note there is no evolution by the Yarkovsky/YORP effect in this simulation. The dotted vertical line denotes the 
value CTe = 0.046 of the observed Hilda family. 



5.2 A sensitivity to the impact-related parameters 

Another important test was devoted to the impact pa- 
rameters, which were varied in a relatively large steps: 
ei = (0.12,0.15), ii = (6.8°, 8.8°), /imp = (45°, 90°, 135°), 
Wimp = (60°, 90°), a = (2.25,4.25), «max = (200, 400) m/s, 
RPB = (83.5, 103.5) km, ppb = (1000, 2000) kg/m^ Note 
that the selection of impact parameters is rather extreme, 
so that we do not expect they may ever be out of these 
bounds. The total number of simulations is 384. The migra- 
tion parameters were fixed (they correspond to one success- 
ful migration scenario): ay, = 5.2806 AU, asi = 8.8250 AU, 
eji = 0.065, esi = 0.08, rmig = 3 Myr, edampj = 10"", 

GdampS — 10 

This time, we decided to use a simple quantity to discuss 



the results, namely the eccentricity dispersion Ue of the syn- 
thetic family at the end of the giant-planet migration. The 
most frequent values of the dispersion are (Je = 0.015 to 0.04 
(see the histograms in Figure [T4|) . Further evolution by the 
Yarkovsky/YORP effect would increase the dispersions up 
to (Te = 0.045 to 0.06, while the observed dispersion of the 
Hilda family is tJe = 0.046. 

We see the histograms look similar for all the impact 
parameters, there is even no apparent correlation between 
them. The explanation for this 'lack of dependence' is that 
the eccentricity distribution is mainly determined by the 
perturbations of the giant planets. A given planetary evolu- 
tion therefore gives a characteristic value of whatever 
the impact parameters are. The dispersion in values 
is due to the fact that the planetary evolutions that we 
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Figure 15. An example of the orbital evolution of Jupiter and 
Saturn with a rare temporary capture in the mutual 3:7 resonance 
(bottom panel). This sort of evolution leads to a large spread of 
pseudo-proper eccentricities of the synthetic "Hilda" family by 
the end of the migration (top panel). 

computed change widely from one simulation to another. 
Though planet migration was prescribed analytically, there 
are mutual interactions of planets and random captures in 
resonances (or jumps across resonances) which may affect 
the eccentricity distribution of the synthetic "Hilda" family. 
An extreme case is shown in Figure 1151 In this particular 
simulation, Jupiter and Saturn were captured in the mutual 
3:7 resonance for 0.5 Myr which resulted in a large eccen- 
tricity dispersion (Te = 0.044 of the synthetic family. Our 
conclusion is that the impact parameters are less important 
than the parameters related to migration. 

5.3 Matching results together 

Even though we do not perform a joint integration which in- 
cludes both the planetary migration and Yarkovsky/YORP 
effect, we try to match the previous results from Sections l5.ll 
and 13.51 together. We do it by using a straightforward 
Monte-Carlo approach: (i) we take the pseudo-proper ec- 
centricities Gmig of bodies at the end of planetary migration 
from Section l5.ll (ii) we compute total Yarkovsky/YORP 
drifts Acye in eccentricity from Section [3.51 (iii) we assign 
every body a drift randomly (cflnai = emig + Aeye) and this 
way we construct an evolved synthetic familylf] Finally, we 
compare the synthetic family to the observed Hilda family 
by computing a Kolmogorov-Smirnov test for A'^(<efinai) and 
A''(<e)obs distributions. 

To avoid problems with low number of bodies (60 in 
case of planetary migration), we perform the above pro- 
cedure 100 times, always with a different random seed for 
the assignment of the Acye. We then take a median of the 
100 KS statistics as a result for one particular run. The 
resulting histograms of the median Dks for various initial 
conditions are shown in Figure [161 

We confirm the conclusions from Section 15.11 — those 
migration scenarios that preserve the largest number of fam- 
ily members (i.e., high A^'ieft) are the same, for which we can 

^ Note that gravitational perturbations, caused by planetary mi- 
gration, are independent of size (mass), so a large body may be 
easily found at the outskirts of the family. This is another reason 
for the random assignment of Yarkovsky/YORP drifts. 



find a good fit of eccentricity distribution (low Dks)- More- 
over, it seems we can exclude also the timescale of migration 
Tmig = 30 Myr since the total number of successful simula- 
tions is significantly smaller in this case. 



6 CONCLUSIONS 

Results of this paper can be summarised as follows: 

(i) The Hilda family evolves mainly due to the 
Yarkovsky/YORP effect and the observed large spread of 
eccentricities indicates the age 4^^ Gyr. 

(ii) The collisional disruption of a D ~ 200 km parent 
body is unlikely in the current environment. Instead, it 
rather occurred during the Late Heavy Bombardment when 
collisions with comets dominated and were up to 50 times 
more frequent. Another possible source of projectiles is the 
population of D-type asteroids captured in the J3/2 reso- 
nance (Levison et al. 2009). 

(iii) In case the Hilda family was created during giant- 
planet migration, which seems to us likely, the major per- 
turbations of the family were due to secondary resonances 
between libration frequency and the frequency of Jupiter- 
Saturn 1:2 critical argument. 

(iv) On the basis of our simulations, we argue the mi- 
gration was relatively fast (with time scale Tmig — 0.3 Myr 
to 3 Myr) and Jupiter and Saturn were relatively closer to 
the current configuration (with period ratio Ps/Pj J? 2.13 
or more) at the moment when the "Hilda" family was cre- 
ated, otherwise the family would be 'destroyed' by migra- 
tion. Slower migration time scales are only allowed for larger 
values of Ps/Pj ratios. 

The Hilda family thus proved to be one of the oldest families 
in the main asteroid belt. 

There are emerging indications that orbital evolution 
of planets was rather violent and close encounters between 
planets were present (Nesvorny et al. 2007, Brasser et al. 
2009). This might be still consistent with our model of the 
Hilda family, but of course we have to assume the family 
formed after severe perturbations in the J3/2 region ended. 
A more complicated migration scenario like that of 'jumping 
Jupiter' (Morbidelli et al. 2010) even seems favourable in our 
case because Jupiter and Saturn very quickly reach a high 
period ratio (Ps/Pj ^ 2.3, i.e. the planets are quite close 
to their current orbits). Then, the perturbations acting on 
the J3/2 region are already small and the fiux of impactors 
becomes high just after the jump. The Hilda family thus 
might have formed exactly during this brief period of time. 

Regarding future improvements of our model, knowl- 
edge of geometric albedos for a large number of small aster- 
oids may significantly help and decrease uncertainties. The 
WISE infrared mission seems to be capable to obtain this 
data in near future. 
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Figure 16. The number of simulations N versus the Kolmogorov— Smirnov distance Dks between the synthetic and the observed Hilda 
family. The simulation differ by the time scale of migration r^ig and the initial conditions for Jupiter and Saturn (aji, a-si)- We only 
plot successful migration scenarios with Aiipi^^ptg ^ 2000 m/s and the number of bodies left A^i^ft > Afinitiai/2. The dotted vertical line 
denotes the distance Dks foi' which the probability p(>£'ks) that the two eccentricity distributions differ by this amount equals to 0.01. 
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